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As was demonstrated in earlier studies, turbulence can result in a negative contribution to the effective mean magnetic 
pressure, which, in turn, can cause a large-scale instability. In this study, hydromagnetic mean-field modelling is performed 
for an isothermally stratified layer in the presence of a horizontal magnetic field. The negative effective magnetic pressure 
instability (NEMPI) is comprehensively investigated. It is shown that, if the effect of turbulence on the mean magnetic 
tension force vanishes, which is consistent with results from direct numerical simulations of forced turbulence, the fastest 
growing eigenmodes of NEMPI are two-dimensional. The growth rate is found to depend on a parameter /?* characterizing 
the turbulent contribution of the effective mean magnetic pressure for moderately strong mean magnetic fields. A fit 
formula is proposed that gives the growth rate as a function of turbulent kinematic viscosity, turbulent magnetic diffusivity, 
the density scale height, and the parameter P^ . The strength of the imposed magnetic field does not explicitly enter provided 
the location of the vertical boundaries are chosen such that the maximum of the eigenmode of NEMPI fits into the domain. 
The formation of sunspots and solar active regions is discussed as possible applications of NEMPI. 
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1 Introduction 

The concept of turbulent viscosity is often used in astro- 
physical and other applications in recognition of the fact 
that the microscopic viscosity is far too small to be relevant 
on the length scales under consideration. Turbulent viscos- 
ity is the simplest parameterization of the Reynolds stress 
tensor, UiUj, where u = U — U is the velocity fluctuation 
about a suitably defined average, denoted here by an over- 
bar. Turbulent viscosity is by far not the only contribution 
to the Reynolds stress tensor. In addition to hydrodynamic 
contributions such as the A effect (Riidiger 1980, 1989), 
which is relevant to explaining stellar differential rotation 
(Riidiger & Hollerbach 2004), and the anisotropic kinetic 
alpha effect (Frisch et al. 1987), which provides an impor- 
tant test case in mean-field hydrodynamics (Brandenburg & 
von Rekowski 2001 ; Courvoisier et al. 2010), there are mag- 
netic contributions as well. One can think of them as a mag- 
netic feedback on the hydrodynamic stress tensor (Radler 
1974; Riidiger 1974) or, especially when magnetic fluctua- 
tions are also considered, as a mean-field contribution to the 
turbulent Lorentz force. 

Work by Roberts & Soward (1975) and Riidiger et al. 
(1986, 2012) using a quasi-linear approach suggests that 
the total magnetic tension force (which includes the effects 
of fluctuations) is reduced in the presence of a mean mag- 
netic field and might formally even change sign for larger 
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magnetic Reynolds numbers, but this would then be be- 
yond the validity of their approximation. Using the spec- 
tral r approach, Kleeorin et al. (1989, 1990) do indeed find 
this reversal of the sign of the total magnetic tension force. 
In addition, they find a reversal of the sign of the effective 
magnetic pressure term; see also Kleeorin & Rogachevskii 
(1994) and Kleeorin et al. (1993, 1996). Rogachevskii & 
Kleeorin (2007) argue that, in a stratified medium, this can 
lead to the formation of large-scale magnetic flux structures 
and perhaps even sunspots - or at least active regions. 

Recently, direct numerical simulations (DNS) of both 
unstratified and stratified forced turbulence (Brandenburg et 
al. 2010, 2012; hereafter referred to as BKR and BKKR, 
respectively) have substantiated this idea and have demon- 
strated that the effective magnetic pressure can indeed 
change sign. Similar results have now also been obtained 
for turbulent convection (Kapyla et al. 2012). In addition, 
these papers give results of mean-field calculations illus- 
trating that there is a negative effective magnetic pressure 
instability (hereafter referred to as NEMPI) when there is 
sufficient density stratification. 

NEMPI is a convective type instability related to the in- 
terchange instability in plasmas (Tserkovnikov 1960; New- 
comb 1961; Priest 1982) and the magnetic buoyancy insta- 
bility in the astrophysical context (Parker 1966). The free 
energy in interchange and magnetic buoyancy instabilities 
is drawn from the gravitational field, while in NEMPI it is 
provided by the small-scale turbulence. 

The mechanism of NEMPI works even under isother- 
mal conditions when entropy evolution is ignored and an 



®W]LEY 

jl-nterScience* 



© 2012 WILEY-VCH Verkg GmbH& Co. KGaA, Weinheim 



96 



K. Kernel et al.: Negative effective magnetic pressure instability 



isothermal equation of state is used. This has been shown 
using corresponding mean-field calculations (BKKR). With 
this reduction to the most elementary aspects of the insta- 
bility, it has recently been possible to verify the existence 
of NEMPI also in DNS (Brandenburg et al. 20 11, hereafter 
referred to as BKKMR). This has been a major step for- 
ward, because now there is no doubt that one is pursuing 
a real effect and not just one that works only in the world 
of mean-field models. Essential to the paper of BKKMR 
has been a finding from an earlier version of the present 
one that only two-dimensional mean-field structures are ex- 
cited. This property allowed meaningful averaging along the 
direction of the imposed field, making the identification of 
flux concentrations thus much clearer The absence of three- 
dimensional mean-field structures was surprising because 
three-dimensional mean-field calculations have shown that 
the mean magnetic field develops structures along the di- 
rection of the imposed field (BKR). However, while the 
mean-field calculations have illustrated the nature of the in- 
stability, no systematic survey of solutions has yet been at- 
tempted. The purpose of this paper is therefore to clarify 
some still puzzling aspects concerning NEMPI. Note also 
that the large-scale flux concentrations observed in DNS 
of BKKMR have an amplitude of only 15 % of the local 
equipartition field. This implies that the flux concentrations 
we observe in DNS are often not strong enough to be no- 
ticeable without averaging. 

In addition to the structures found in BKKMR, other 
types of structures have recently been reported in Large- 
Eddy Simulations (LES), which might also be an indication 
of NEMPI. We have here in mind the radiation magneto- 
convection simulations of Kitiashvili et al. (2010), in which 
one sees the formation of whirlpool-like magnetic struc- 
tures. Relevant to NEMPI is also the work of Tao et al. 
(1998), who considered magneto-convection in the optically 
thick approximation and find a horizontal segregation into 
magnetized and non-magnetized regions. The size of the in- 
dividual regions is such that they encompass several turbu- 
lent eddies. This phenomenon might therefore well be as- 
sociated with an effect that could also be modelled in terms 
of mean-field theory. However, before we can make such an 
association, we need to find out more about the properties 
of NEMPI. In particular, we need to know what is the op- 
timal magnetic field strength, what are the requirements or 
restrictions on the turbulent velocity, and, finally, how much 
density stratification is needed to make NEMPI work. 

To connect the aforementioned requirements to DNS, 
we need to have a meaningful parameterization of the tur- 
bulence effects. The work done so far has been focussing 
on measuring a reduction of the turbulent pressure and ef- 
fective mean magnetic pressure as a function of the local 
mean magnetic field strength. The shape of the resulting de- 
pendence of the effective mean magnetic pressure on the 
mean magnetic field has been matched to a specific fit for- 
mula that can be characterized by two fit parameters that, in 
turn, can be linked to the minimum effective mean magnetic 



pressure and critical field strength above which the effect is 
suppressed. However, there have been indications that this 
parameterization is not unique and that different combina- 
tions of the two fit parameters can result in similar values of 
minimum effective pressure and the critical field strength. 
The question therefore arises whether this apparent degen- 
eracy affects the properties of NEMPI. 

We mentioned already the fact that NEMPI is capable 
of exciting three-dimensional structures that show varia- 
tion along the direction of the mean magnetic field. This 
would give rise to the worry that the two-dimensional re- 
sults presented so far may not reflect the properties of the 
fastest growing mode and may therefore not be relevant to 
describing NEMPI. However, as will be discussed in this 
paper, this is not the case, because the degree to which 
three-dimensional modes are excited depends on the sign 
and magnitude of one of the turbulence parameters, namely 
the term characterizing turbulence effects on the magnetic 
tension force, and that simulations indicate that this sign is 
not favorable for exciting three-dimensional modes (BKKR, 
Kapyla et al. 2012). Before we begin addressing the vari- 
ous points, we discuss first the mean-field model and basic 
setup. 

2 Mean-field model 

In view of further verifications of NEMPI with DNS, it is 
necessary to be able to reduce the essential physics to a 
minimum. We will therefore not make any attempt to con- 
sider other aspects that would make the model more realistic 
with respect to the Sun. Given that NEMPI works even un- 
der isothermal conditions (BKKR), we adopt an isothermal 
equation of state where the mean pressure p is linear in the 
mean density p, with p — 'pc^ and Cg being the constant 
isothermal sound speed. We solve the evolution equations 
for mean velocity U, mean density p, and mean vector po- 
tential A, in the form 

ffFj 

— = -U-VU-clVlnp + g + J^M+J^K, (1) 

^ = -U-Vp-pV-U, 
dA 



(2) 



dt 



UxB-ir^t+v)J, 



(3) 



where !Fm is given by 

p^M = -iV[(l - qp)W] +B-V[{1- qs)B\, (4) 

and 

^K = (t^t + !^)(v2[7+iVV-C7 + 2SVlnp) (5) 

is the total (turbulent plus microscopic) viscous force. Here, 
Sij ~ ■^{Ui,j + Uj^i) — |<5ij V • U is the traceless rate of 
strain tensor of the mean flow. As in earlier work (BKR, 
BKKR), we approximate ^p and q^ by simple profiles that 
are only functions of the ratio /3 = |i?|/i3oq. However, in 
the earlier work this functional form was described by 

qaifi) - ?.o[l - (2/^) arctan(/3V/3')], (6) 
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where B^q is the equipartition field strengths and a stands 
for subscripts p and s, respectively. We refer to this as the 
arctan fit. In the present paper we use an algebraic fit of the 
form 



qa{(3) = q.o[l - (2/^) arctan(/?V/3')], 



qa{(3) = 



1 + /97/3; 



2 ■ 



(7) 
(8) 



The functions qp and Qs quantify the impact of the mean 
magnetic field on the effective pressure and tension forces, 
respectively. 

As initial condition, we assume a hydrostatic stratifica- 
tion with p(z) = po exp(— z/i/p), where Hp ~ c^/g is the 
scale height in our domain of size Lx^^LyXLz', the exact 
dimensions vary between 4 and 10 density scale heights in 
each direction. We add a small perturbation to the velocity 
field. We allow for the presence of an imposed field in the y 
direction, Bq = (0, Bq, 0). The total field is then written as 

B = Bo + V xA, (9) 

so the departure from the imposed field is expressed in terms 
of the mean magnetic vector potential A. Furthermore, we 

assume 



B. 



cqiz) = Boqoexp(-z/2iJp), 



(10) 

with a normalization coefficient i?cqo- This formula is com- 
patible with Bcq = p^^'^Mrms in BKKMR, where the turbu- 
lent rms velocity, Umis, was approximately constant. 

On the upper and lower boundaries we adopt stress-free 
boundary conditions for velocity, i.e. Ux,z — Uy^z = Uz = 
0, and a perfect conductor boundary condition for the mag- 



netic field, i.e. Ax = Ay = Az 



0. Here, commas denote 



partial differentiation. No boundary condition for the den- 
sity is required. All computations have been carried out with 
the Pencil Cod^ 

Our model is characterized by the following set of input 
parameters. There are three parameters characterizing the 
hydrostatic equilibrium stratification, namely g, Cs and po. 
The remaining parameters are the normalized imposed field 
strength, Bo/Bcqo, turbulent viscosity and magnetic diffu- 
sivity, as well as the parameters Qo^q and /S^- 

3 Results 

3.1 Two- and three-dimensional solutions 

Earlier work has suggested that the eigenmodes of NEMPI 
can be three-dimensional (BKR). This could render two- 
dimensional calculations inadequate if the first excited 
mode were indeed three-dimensional. However, it turns out 
that the wavelength of the eigenmode in the direction of 
the field increases as Qs decreases. In BKR, where three- 
dimensional (y-dependent) solutions to NEMPI were first 
reported, Qs was chosen to be around 10, and the fastest 
growing mode was indeed three-dimensional. In Fig.[T|we 
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Fig. 2 Dependence of ky on qso (upper panel), together 
with the corresponding growth rate A (lower panel). 

Table 1 Comparison of normalized growth rates, 

XHp/rjt, for different values of qs, for a three (3D) and two- 
dimensional (2D) simulation (Ly -^ cx)). 
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show that the effective wavenumber of the variation of the 
field in the y direction decreases with decreasing values of 
qs- This is shown quantitatively in Fig. |2l where we plot 
the dependence of the typical value of the field-aligned 
wavenumber, ky, on the value of QsO- Here, ky is evaluated 
in a layer near the surface. 

We find that the typical value of ky grows with increas- 
ing values of ^so- In addition, we find that the growth rate 
of the instability. A, increases with qso approximately Un- 
early once qso exceeds a value of around two. The fact that 
ky ^ as qso — J^ is significant, because BKKR and also 
Kapyla et al. (2012) found from simulations that qso ~ 0. 
In that case, the characteristic length scale along the direc- 
tion of the field becomes infinite and the calculation essen- 
tially two-dimensional. Conversely, when studying NEMPI 
in two dimensions, changing the value of q^Q has no ef- 
fect on structure formation and the growth rate; see Table [T] 
However, it is now clear that this is an artifact of restricting 
the solutions to be two-dimensional. 



3.2 Approximate degeneracy in the Qp fit formula 

We mentioned in the introduction that recent attempts to de- 
termine qpQ from simulations faced the difficulty that the fit 
formula (|7j possesses an approximate degeneracy in that we 
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g=o-2o 




Fig. 1 (online colour at: www.an-journal.org) Visualization of By at the periphery of the computational domain near the 
end of the kinematic growth phase. Note the change of the field pattern with increasing values of Qsq — 2, 5, 10, and 20 
(from left to right). 
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Fig. 3 Top: comparison of original fit curves with differ- 
ent (7po values (10, 20, 50, 100, 300) and a fixed minimum, 
showing similar zero crossing values for a large qpo range. 
Bottom: same as above, but now using an algebraic fit, giv- 
ing larger spacing between the curves. 



It is therefore better to measure the normalized minimum ef- 
fective magnetic pressure, "Pmin — \ min[(l — gp)/3^], and 
in- For the algebraic fit we then obtain the fit 



its position, fi-a 
parameters 

/^P = ^min / V 



-2P„ 



/?* = /3p + V-SK, 



(12) 



where we have introduced the parameter /3j — gpo/3p in a 
modified representation 

PI 



9p(/3) 



PI 



./32 



(13) 



which is preferable over Eq. (O in circumstances where 
01 = q-pQ0i is approximately constant. This appears to be 
the case in recent DNS (BKKR, Kemel et al. 2012), where 
/?v, « 0.2 and 0.3 in the absence and presence of small-scale 
dynamo action, respectively. 

We have computed mean-field models for different com- 
binations of parameters using the algebraic fit. We find that 
the resulting growth rate A depends on the functional form 
of T^cff (/?) near /3 = /^min, which manifests itself in a de- 
pendence on both (7po and /3p; see Fig. |4] The lower panel 
of this figure suggests that the dependence of the growth 
rate on both parameters can be collapsed onto a single de- 
pendence on /?*. This underlines the usefulness of Eq. ( flj] ) 
as a fit formula. As argued above, this dependence is best 
constrained by the fit parameters /Smin and 'Pniin- 



can obtain a similarly looking dependence of the effective 
mean magnetic pressure, 

7'cff(/3) = i[l--Zp(/3)]/32, (11) 

over a wide range of values of ^po by adjusting the value 
of /3p correspondingly. This can be seen in Fig. [S] where 
we show Vcs{P) using either the arctan fit (upper panel) or 
the algebraic fit (lower panel) for parameters that result in 
the same value of "Pniin for a range of values of q-^o. Note 
that the position where VcS becomes positive, i.e. the crit- 
ical value defined by "Pcffl/^ciit) = 0, is rather similar in 
all cases. This approximate degeneracy is particularly obvi- 
ous for the arctan fit, and less so for the algebraic fit. How- 
ever, in both cases the form of VcS near /3 — > changes 
significantly. Therefore, the approximate degeneracy would 
be lifted if one could determine q-pQ from the behavior of q-p 
near (3 — 0. However, near /? = the DNS have large errors. 



3.3 Onset of NEMPI 

With a given prescription of qp{(3), assuming here qs — 0, 
we can now compute two-dimensional mean-field models. 
Our goal is to obtain a simple formula that can tell us how 
large the growth rate of the instability is, and what the crit- 
ical condition for the onset of the instability is. Not much 
is known about the linear stability properties of NEMPI, so 
we have to rely on numerical determinations of the growth 
rates for different wavelengths for different parameters to 
obtain an approximate representation of the dispersion re- 
lation. Earlier work of Kemel et al. (2011) has suggested a 
relation of the form 

A = ^{g/clqpo, Pp, ...) ~ vtkl - 77tA;2, (14) 

where fc^ and A;,, are effective wavenumbers quantifying the 
effects of turbulent viscosity and turbulent magnetic diffu- 
sivity, $ is a function of the inverse scale height, Hp^ = 



© 2012 WILEY- VCH Verlag GmbH & Co. KGaA, Weinheim 



www.an-joumal.org 



Astron. Nachr. / AN (2012) 



99 



100 




100 r 



0.01 



0.10 



100 



10 




0.1 



1.0 



Fig. 4 Dependence of the normalized growth rate A/Aq, 
with Ao = (I't + 7?t)/^p, on qpo for /3p = 0.05 (up- 
per panel), on /3p for qpo = 40 (middle panel), and on 

1/2 

13-1, = gpQ /3p (lower panel) for /3p = 0.05 (red), gpo = 40 
(circles), and /3p — 0.1, Qpo = 5 (blue). Solid lines represent 
approximatefits given by 9(gpo — l)°'^, (/3p/0.0015)^-^, and 
(/3,,/0.02)3/2, respectively. 
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Fig. 5 

i 
fed ~ ^.h/Hp, where Hp 



ately in terms of g and Cg. The negative slope gives fc^ with 



5/Cg, and other parameters describing the functional form 
of(?p. 

We now need to determine the various unknowns. We 
begin by determining ki, and fc^ by varying either only vt 
or only r]t at a time. It turns out that k^ = fc,, = fed suf- 
fices. In this way we obtain a linear fit for the growth rate 
A = const — (vt + ?7t)fcd' gi'^iiig ^^ ^d ^^ '■^^ slope of this 
graph; see Fig.|5] We find fed ~ 2.5/ Hp, superseding earlier 
results by Kemel et al. (201 1) for a different Bcq{z) profile. 
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Fig. 6 Dependence of $ = A+(t't+r?t)fcd o" (feiJp)^^/"' 
(upper panel), A on Bo/Beqo (middle panel), as well as 
zb/ Hp and z-p/Hp on Bo/Bcqo (lower panel). Here, /3^o = 
0.0083 is a fit parameter and Ad = {t^t + Vt)k'^ is used for 
normalization. 



Accepting now the fit parameter fed as measured, we can 
proceed to determining the dependence of <I> on Hp (top 
panel of Fig.|6]l. A convenient non-dimensional quantity is 
k{Hp, where fcf is the wavenumber of the energy-carrying 
eddies of the turbulence, which is related to r/t — Umis/Sfcf 



with Mr 



B, 



cqO 



/Po 



1/2 



Note that $ ex (fcfi/p)"^/'*. Com- 
bining this with the /3j scaling of Fig. 21 we suggest 



A 



(/3.//3.C 



)3/2(fcfi/^)-3/4 



ii't+Vt)kl 



(15) 



with /3-to w 0.008 being yet another fit parameter Inter- 
estingly, A is independent of the imposed field strength, 
Bo/B^qp, provided the bulk of the eigenmode (z = zb) 
fits well within the domain (middle panel of Fig. |6]). This 
has here been achieved by adjusting the positions of the 
boundaries, ztop and Zbot- Indeed, as Bo/B^qo is increased, 
Zb is found to decrease approximately like ZB/Hp « 
— 21ni3o/-BcqO + const. It turns out that zb is about 2-3 
scale heights below the location z-p where VeS^z) attains 
its minimum value. Contrarily, when z-p > Zbot > zb or 
when Ztop < zb NEMPI will depend on respectively /3bot 
or /3t,op, when zp < Zbot there is no instability. 

Some comments about the horizontal dimensions are 
in order In all cases with gso = 0, we find that in three- 
dimensional calculations with finite y extent, the value of 
Ly does not affect the growth rates. On the other hand, dou- 
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bling the x extent yields two pairs of rolls. This has also 
been confirmed for DNS; see Kernel et al. (2012). 

4 Conclusions 

The present work has clarified a number of puzzling as- 
pects of NEMPI. Firstly, it is now clear that we can proceed 
with two-dimensional mean-field simulations as long as we 
know that gso = (or negative). However, this may not 
always be the case. The fact that three-dimensional struc- 
tures can emerge from NEMPI was initially thought to be 
an interesting aspect, because it could readily explain the 
formation of bipolar regions (BKR). However, given that 
simulations now indicate that qs ~ (or perhaps even neg- 
ative), this proposal would thus not be an option, unless 
some other as yet unexplored effect begins to play a role. 
In principle, all turbulent transport processes are nonlocal 
and must be described by a convolution with the mean field 
rather than a multiplication (Brandenburg et al. 2008). In 
Fourier space, the convolution corresponds to a multiplica- 
tion with a wavenumber-dependent turbulent transport co- 
efficient. Thus, the idea of explaining bipolar regions would 
again become viable if this effect only existed at small and 
intermediate length scales. Clarifying this would be a task 
for future simulations, because none of the currently avail- 
able techniques are yet equipped to address this possibility. 

Next, we have seen that the degeneracy in the fit for- 
mula used for qp{(3) and Pcff (/3) is significant in that differ- 
ent combinations of qpo and /3p result in similar values of 
min(Pcff ) and /3crit, but the growth rates can still be quite 
different. This means that it is not sufficient to measure only 
min(7'cff) and /3crit- Instead, to characterize the functional 
form of ■Pcff (/3) more accurately, we need some other char- 
acteristics to represent the dependence of this function near 
13 = 0. One such possibility is to use the field strength (3nun 
for which the minimum of the effective magnetic pressure 
is reached. 

Knowing the value of /?min has particular relevance 
in determining the height where NEMPI occurs. For a 
given value of the imposed field strength Bq, the condition 
Bo/Bcq{z-p) = /3min determines the height z-p, where the 
effective magnetic pressure attains a minimum, and thus the 
height zb, which tends to be 2-3 scale heights below Zmin; 
see the lower panel of Fig. |6] Therefore, the value of Bq 
does not directly affect the growth rate of NEMPI. 

Finally, we have tried to establish an approximate dis- 
persion relation to estimate the growth rate of NEMPI as 
a function of turbulent viscosity, turbulent magnetic diffu- 
sivity, mean field strength, and the strength of stratification. 
This formula may serve as a first orientation and can hope- 
fully be improved further with future simulations. This for- 
mula can also be useful in connection with analytic esti- 
mates concerning the regimes when NEMPI is expected in 
DNS or under other more realistic circumstances. 
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